Density-matrix functionals for pairing in mesoscopic superconductors 
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A functional theory based on single-particle occupation numbers is developed for pairing. This 
functional, that generalizes the BCS approach, directly incorporates corrections due to particle 
number conservation. The functional is benchmarked with the pairing Hamiltonian and reproduces 
perfectly the energy for any particle number and coupling. 
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While the Bardeen-Cooper-Schrieffer (BCS) micro- 
scopic theory [l[ provides a suitable description of super- 
conductivity in the macroscopic limit, it is not accurate 
enough for mesoscopic systems, such as nuclei, atomic 
clusters, quantum dots, fullerenes, nanotubes, or ultra- 
small metallic grains ^2-^. Standard BCS approach to 
superconductivity has some dravifbacks: (i) non negligible 
corrections due to finite size effects are necessary in meso- 
scopic systems; (ii) in condensed- matter and/or nuclear 
physics, BCS theory starting from the bare many-body 
interaction cannot be considered either as a numerically 
tractable approach nor as a predictive theory. To over- 
come difficulty (ii) , specific functional theories based not 
only on the local density pir) but also on the anomalous 
density K{r,r') are used [6, 7]. Guided by the Hamilto- 
nian case, the energy is decomposed as: 



(1) 



where £ [p] is often taken as the energy density functional 
without pairing, while £pair[K, k*] is the extra contribu- 
tion due to pairing. Eq. ([!]) turns out to be very accu- 
rate to deal with (ii) for instance in nuclear physics. In 
that case, systems with 10-200 constituents are consid- 
ered and additional finite size corrections are necessary 
0, 0, @]- However, recent studies have shown that tech- 
niques generally used to restore good particle number 
should be handled with care when combined with den- 
sity functional theories [ol-lisj. In particular, unless new 
methods able to properly treat finite size effects and more 
generally configuration mixing within functional theory, 
most of the functional designed during the last 30 years 
have to be revisited [l^. These difficulties question the 
possibility to use symmetry breaking within a functional 
theory. 

The goal of the present work is to provide a new theo- 
retical framework for pairing in finite systems that avoids 
difficulties recently encountered in functional theories 
and that can be applied easily be implemented in cur- 
rent functionals. Here, we follow the idea of Gilbert [isl- 
17[ and seek for a functional of occupation probabilities 
Hi and natural orbitals fi, i.e. £ = £[ni,tpi]. Note in 
passing, that ([T]) already enters into the class of Gilbert 
functionals. Indeed, the pairing energy is generally writ- 



ten as £pair[K, K*] = J2i]K]ki'^ij'^ki where v'^'^ denotes 
the effective two-body kernels in the particle-particle and 
hole-hole channels. In practice, both p and k are com- 
puted using a quasi-particle (QP) vacuum trial state 
that can be seen as a generalization of the Kohn-Sham 
Slater determinant. If the energy is minimized in the 
canonical basis, the trial state, denoted by |(/)), takes 

a BCS like form \(j)[ni, ipi]) = Hi (l + Xiolal) |0) with 



Xi — y^ni/{l — rii) and where |0) corresponds to the 
particle vacuum while {a|,at} are associated to dou- 
bly degenerated canonical states {Lpi,Lpj\ with occupa- 
tion probabilities n^. In the canonical basis, k becomes 
block diagonal with k^j — -^Z 71^(1 — ni) and finally leads 
to an energy functional of the BCS occupation numbers. 
Functional based on BCS suffers for instance from the 
absence of pairing at weak coupling, it also misses part 
of the pairing effects at strong coupling (see for instance 
[l3|)- These defects can be cured by considering many- 
body trial states projected onto good particle number. In 
that case, a reference QP state is first introduced, onto 
which the projection is made. The resulting energy be- 
comes a rather complicated functional of the reference 
state [3]. However, serious difficulties appear when pro- 
jection technique made before or after variation is com- 
bined with functional theories 

Here, we use a completely different strategy and con- 
sider the projected state directly as the trial state. This 
state, called hereafter. Projected BCS state is written in 
its canonical basis as: 



\N) 



1 



(rt)A'l-) 



(2) 



where = Xib\ with b\ — a\a]. N denotes here the 
number of pairs. From this state, we define the occupa- 
tion number and two-body correlation matrix elements 
through: 
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(3) 



Guided by the Hamiltonian framework, we propose to 
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generalize the pairing energy and writes 



(4) 



where w"^" can eventually depend on the density of the 
projected state. By doing so, difficulties observed with 
projection are avoided. The main result of the present 
paper is to show that the correlation can be accurately 
written as a functional of occupation numbers of the pro- 
jected state and account for finite size correction. 

According to the definitions ([3]), both n. 
be written as functionals of the parameter set {xi} 



and Rij can 
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for (iy^j), 



(5) 



while for i — j, Rij — rii. In these expressions, the 
following quantities have been defined (for any K such 
that 1 < K < N and any n with n < H.): 



^3K I 



,(6) 



where means that the summation is made over all 
(jii ■ ■ ■ Jk) different from each other while (ji, • • ■ jx) 7^ 
• ■ ■ I *n) adds the constraint that all j„ are different 
from (zi, • • • , i„). From the above discussion and as can 
be intuitively from the expression of the trial state (dJ, 
the energy can be written as an explicit functional of 
the {xi}. Unfortunately, the complexity of this func- 
tional prohibits its use^. A second difficulty of using 
{xi} as variational parameters is that they are not easily 
connected to quantities like occupation numbers adding 
complexity in physical interpretation. The BCS exam- 
ple illustrates that the energy can directly be written as 
a functional of the {ui} but is insufficient to precisely 
grasp the physics of pairing in finite systems. 

Despite the complex relations between the occupation 
numbers and the {xi} (Eq. (O), it is shown below that 
(i) the {xi} can be accurately replaced by a functional 
of the {n-i}; (ii) using this functional, the energy itself 
becomes a functional of the occupation probability only, 
that provides a very accurate description of the pairing 
Hamiltonian energy. Different relations given below are 
directly derived from recurrence relations existing for the 
jR- |19l42lj after tedious but straightforward calculations. 
We only give here main results are summarized here while 
technical details will be given elsewhere j2l| . 



Eqs. ^ can be rewritten as 

I n -^i -^j / _ \ fry\ 

\Xi\^ + ai \XiY-\Xj\^ 

where Ui = / {N In 

An interesting connection with the BCS theory can be 
made by inverting the first expression in ([7]) as 



1 - n. 



(8) 



Reporting in the second equation of ([7]), leads to 



Ri 



(9) 



Assuming ai — aj — 1, this expression identifies with the 
BCS limit. Therefore, the physics associated with parti- 
cle number restoration is all contained in the set of {ai} 
parameters. The a parameters can be written as a func- 
tion of the occupation numbers by performing a (l/N) 
expansion assuming that the leading order identifies with 
the BCS hmit one obtains^ [2l| : 
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(10) 



which is nothing but a; written as an explicit functional 
of the occupation numbers. By replacing a; given above 
in R and then the expression of i? in the energy 
become an explicit functional of the occupation numbers 
through the sequence 

£[p,R,(pi] -> £[{xi},ipi] ^ £[{ni},(p.i]. (11) 

In practice, the new expression we obtained for £ can di- 
rectly replace the BCS expression in theories either based 
on a Hamiltonian or directly formulated in a density func- 
tional framework. Since the functional is directly formu- 
lated in terms of natural orbitals and occupancies, one 
should add specific constraint during the minimization. 
Below, the quantity: 

£'[{n,},^,] - £[{n,},ip,]^fi(y2n,~N), (12) 



^ Note that recently, it has been shown that such a minimization 
could eventually be performed numerically using specific recur- 
rence relations satisfied by the Ix [19142 ill . 



^ Note that, additional terms tested numerically as negligible and 
appearing at the second order approximation are omitted here. 
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FIG. 1: (Color online) Difference between the Hartree-Fock 
energy £hf and the exact Richardson solution (red solid line) 
of ref. [l^, generally referred as the condensation energy, 
compared to result obtained with successive approximation 
using expression pop for 8 particles. The dashed, dotted, 
dot-dashed lines and filled circles correspond respectively to 
the BCS, first, second and third order correction. In insert, 
a focus on the weak coupling region is shown. Note that the 
PBCS energy (not shown here) cannot be distinguished from 
the exact result. 



pling regime. It could indeed be shown that all terms in 
the expansion up to order 1/7V! contribute equally in the 
Hartree-Fock limit. This directly stems from the inade- 
quacy of the BCS functional in the weak coupling limit. 
For instance, a direct use of perturbation theory leads 
to much better results underlying the role of 2 particles- 
2 holes excitations (see for instance discussion in jl9|). 
This difficulty could only be overcome by considering 
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Va. Va^ 



is directly minimized. 

In order to test the functionals developed here, we con- 
sider an even system of A particles interac ting through 
the "picket" fence Hamiltonian of the form 22 1 



iJ = ^ ej(a|ai -I- ata-) 



g a]ata:;a 
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(13) 



with doubly degenerate equidistant levels and with a level 
spacing Ae. In the following, N = A/2 denotes the num- 
ber of pairs while Vl = AN is the single-particle space 
size. For this Hamiltonian, the energy simply writes 
£(ni, Rij) = 2 EiTii — g Rij- Note that, this func- 
tional interaction corresponds to dH) with v''-''.-. = — 
g. The pairing Hamiltonian can be solved exactly and 
therefore is particularly suitable for benchmarking ap- 
proximation for pairing correlation (sl. [23| . 

An illustration of successive corrections beyond the 
BCS approximation is given in figure [1] Results are ob- 
tained using the functional form of R where ai are re- 
placed by Eq. ^TO\i trunctated at a given order. Then, 
Eq. (jl2p is directly minimized starting from the BCS so- 
lution and making variation of occupation numbers be- 
tween and 1. A quadratic programming (QP) method is 
used for the minimization leading to rapid convergence 
tested up to 400 particles. Systematic improvement is 
observed as higher orders in the correction are included. 
The above functional has however two major drawbacks. 
First, it is rather complicated to use. Second, while few 
terms are necessary to get a perfect result in the strong 
coupling, the convergence is rather slow in the weak cou- 



FIG. 2: (Color online) Comparison between the exact con- 
densation energy (solid line), the BCS (dashed line) and the 
result obtained by minimizing the density matrix functional 
(filled circles) with the approximation (|15l) for various parti- 
cle numbers and coupling strength . Again in all cases, the 
PBCS result is in perfect agreement with the exact one. 

these terms explicitly, which in the form given by ^TU\\ 
becomes impractical for numerical implementation as N 
increases. A simple linear expression, i.e. 



Oil = ao + aiUi, 



(14) 



can however be found by considering all terms in the 
expansion and approximating sums in (jlOp according to: 
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A straightforward calculation then gives: 
1 1 



ai 



^2 



TV 1 - S2 



dai 

ao = 1 - (S2 - S3)-— , (15) 
as2 



where Sp — j^^iiniY ■ Note that, due to the re- 
summation of expansion ([ro]) of all terms up to order 
N , the present approach could not be anymore regarded 
as a 1/N correction. Energies obtained with this approx- 
imation are shown in figure[2]for various particle numbers 
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and couplings. The energy found by minimizing the func- 
tional is overall in very good agreement with the exact 
energy for any particle number and coupling. A careful 
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FIG. 3: (Color online) Evolution of ni(l — n^) as a function of 
Ei/Ae for the exact (solid line), PBCS (open triangles), and 
the new functional (filled circles) for A — 16 and g/Ae — 0.16 
(left) and 0.56 (right). 

analysis shows a slight underestimation of the energy in 
the intermediate coupling regime associated also with a 
slight difference in the occupation numbers (see figure [3|). 
This small discrepancy stems from the linear approxima- 
tion made for the ai. To better account for occupation 
number behavior, quadratic or cubic corrections to (jl4[) 
might eventually be obtained. However, this will add 
complexity to the functional while the energy is already 
rather well reproduced. 

In this letter, a new approach is proposed to account 
for pairing in finite systems using functional theories. To 
escape difficulties recently observed [ol-illl, Il4| . namely 
divergences and jumps, the introduction of an auxiliary 
QP state is avoided and a Projected BCS trial state is 
directly used. In the pairing hamiltonian, a perfect agree- 
ment between the functional result and the exact energy 
is obtained at all coupling strength and particle num- 
ber. The present method can be directly implemented 
on existing functional theories and should provide an ac- 
curate way to treat finite size effects. Note that present 
approach can easily be extended to odd systems [21| and 
might provide a tool of choice to study dynamics and 
thermodynamics of finite systems with pairing. Guided 
by the BCS theory, system at finite temperature can be 
studied minimizing the free energy defined through ^: 



T[{ni}, ipi\ = £'[{ni}, cpi] - TS[ni]. 



(16) 



Note that, the present theory has been recently applied to 
nuclei showing that it solves recent difficulties associated 

^ since all the information on the system is now contained in the 
single-particle components, the entropy simply identifies with 



to broken symmetries 0, [l3| opening new perspectives. 
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